%plot_station_map
network = "HV";
stationlist = ["NPB","SRM","NPT","OBL","WRM","SDH","UWE","UWB",...
    "SBL","KKO","RIMD"];
locations = "*";
channellist = ["HHZ"];

%table names for each channel
channel_station_str = strings(1,length(stationlist)*length(channellist));
ctags = struct(length(stationlist)*length(channellist),0);
for statin = 1:length(stationlist)
    for chanin = 1:length(channellist)
        in = length(channellist)*(statin - 1) + chanin;
        channel_station_str(in) = strcat(networks,'_',...
            stationlist(statin),...
            '_', channellist(chanin));
        ctags(in).network = network;
        ctags(in).station = stationlist(statin);
        ctags(in).channel = channellist(chanin);
    end
end

[ctags,keepind] = local_HV_locations(ctags);

%DEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%import dem
DEM = struct();
[DEM.Z,DEM.x,DEM.y] = import_hawaii_dem();
domain = [257,265,2143.5,2151.5]*1000;
DEM.Z = DEM.Z(:,DEM.x>=domain(1)&DEM.x<=domain(2));
DEM.x = DEM.x(DEM.x>=domain(1)&DEM.x<=domain(2));
DEM.Z = DEM.Z(DEM.y>=domain(3)&DEM.y<=domain(4),:);
DEM.y = DEM.y(DEM.y>=domain(3)&DEM.y<=domain(4));


%downsample
DEM.Z = imgaussfilt(DEM.Z,2);
DEM.downsamp = floor(length(DEM.x)/100);
DEM.x = DEM.x(1:DEM.downsamp:end);
DEM.y = DEM.y(1:DEM.downsamp:end);
DEM.Z = DEM.Z(1:DEM.downsamp:end,1:DEM.downsamp:end);

% disprange = 2*[-max(abs(mogi_HVO_okadastruct.s{RIMDindex}.vert)),max(abs(mogi_HVO_okadastruct.s{RIMDindex}.vert))];
% quivscale = 0.5*range(DEM.x)/1000/...
%     max([max(abs(mogi_HVO_okadastruct.s{RIMDindex}.east)),...
%     max(abs(mogi_HVO_okadastruct.s{RIMDindex}.north))]);
vertcmap = cmocean('balance',64);
% vertcmapedges = linspace(disprange(1),disprange(2),size(vertcmap,1)+1);
% vertcmapedges = [-inf,vertcmapedges(2:end-1),inf];

mu = 30*10^9; %shear modulus
nu = 0.25; %poisson's ratio

epilatlon = [19.405402 -155.281041 100.0];
zone = utmzone(epilatlon(1),epilatlon(2));
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
[epiloc_east, epiloc_north, epiloc_vert] = mfwdtran(utmstruct,epilatlon(1),epilatlon(2),epilatlon(3));